Large-scale second RPA calculations with finite-range interactions 
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Second RPA (SRPA) calculations of nuclear response are performed and analyzed. Unlike in most 
other SRPA applications, the ground state, approximated by the Hartree-Fock (HF) ground state, 
and the residual couplings are described by the same Hamiltonian and no arbitrary truncations 
are imposed on the model space. Finite-range interactions are used and thus divergence problems 
are not present. We employ a realistic interaction, derived from the Argonne V18 potenial using 
the unitary correlation operator method (UCOM), as well as the simple Brink-Boeker interaction. 
Representative results are discussed, mainly on giant resonances and low-lying collective states. The 
focus of the present work is not on the comparison with data, but rather on technical and physical 
aspects of the method. We present how the large-scale eigenvalue problem that SRPA entails can 
be treated, and demonstrate how the method operates in producing self-energy corrections and 
fragmentation. The so-called diagonal approximation is conditionally validated. Stability problems 
are traced back to missing ground-state correlations. 
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I. INTRODUCTION 

Extended RPA theories such as second RPA 
(SRPA) which go beyond first-order RPA, are of- 
ten used to describe the strength, decay width and fine 
structure of nuclear giant resonances (GRs) and other 
collective low-lying states. They also appear particularly 
useful in the case of unitarily transformed realistic in- 
teractions, which have not been calibrated for first-order 
RPA methods, but rather rely on extended model spaces 
to provide converged results. Such interactions also make 
"self-consistent" extended-RPA applications possible, in 
the sense that the ground state and the residual cou- 
plings can be described by the same Hamiltonian. First 
applications using a renormalized Argonne V18 potential, 
derived with the unitary correlation operator method 
(UCOM) were presented in Ref. [|. 

Self-consistent extended-RPA applications in large 
spaces without arbitrary truncations can be envisioned, 
in principle, with any properly constructed finite-range 
effective interaction. However, they are hardly ever 
performed for various technical and conceptual rea- 
sons. Phenomenological effective interactions are fitted 
to sets of experimental data using mostly Hartree-Fock(- 
Bogoliubov) and selected (quasi-particle) RPA results. 
Their range of applicability is then restricted to the se- 
lected observables and many-body methods. Zero-range 
effective interactions, which greatly simplify numerical 
applications, cannot be employed in second-order RPA 
methods, because they do not provide a natural cutoff 
in momentum space. Consequently, no effective interac- 
tions have been fitted to calculations beyond first-order 
RPA and consistency in the treatment of the ground and 
excited states is ordinarily abandoned in practical appli- 



cations of such methods. 

It is the purpose of the present work to perform and 
analyze large-scale (i.e., without arbitrary truncations), 
"self-consistent" (i.e., with a single interaction as the sole 
input) SRPA calculations. We employ mostly the VijCOM 
interaction used in Ref. [J]. It is derived from the Ar- 
gonne V18 potential by means of a unitary transforma- 
tion, which renormalizes it, while preserving the phase 
shifts and retaining the complex structure of the realistic 
interaction. We also use the Brink-Boeker interaction, 
Vbb [5|> which is a simple, central, phenomenological ef- 
fective interaction. No explicit three-body force is used 
at this point. We will not focus on producing realistic 
results and comparing them with data, as was done in 
Ref. [1], but rather on technical and physical aspects of 
the method. We present how the large-scale eigenvalue 
problems that SRPA involves can be treated, demon- 
strate how the method operates in producing self-energy 
corrections and fragmentation, and discuss consistency 
and stability problems. 

In the next section we present the SRPA formalism and 
in Sec. IHII the methods we have used to solve the SRPA 
eigenvalue problem. In Sec. IIVI we discuss our results 
with the help of illustrative examples. We conclude in 
Sec. M 



II. SECOND RPA FORMALISM 

In the following we assume a nuclear Hamiltonian con- 
sisting, in general, of a one-body part and a two-body 
part, 
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Three-body terms are not included. If the total Hamil- 
tonian is considered, then 



A 



2m ^— ' 



is the total kinetic energy, while H 2 = V contains the 
interactions of particle pairs. If the intrinsic Hamiltonian 
is considered, then H 1 = and H 2 = T int + V includes 
the intrinsic kinetic energy of the system, 
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We will employ the SRPA as it was formulated in 
Ref. @ in analogy to RPA. The derivation is based on the 
equations-of-motion method and relies on a quasi-boson 
approximation. We will consider closed- (sub)shell spher- 
ical nuclei and their excited states of definite angular 
momentum and parity J 77 . Excited states are expanded 
in the space of particle-hole (ph) and two-particle-two- 
hole (2p2h) configurations. The symbol p (or h) will 
represent all the quantum numbers of a particle (hole) 
state except the magnetic quantum number m p (rrih), 
i.e., the set of quantum numbers {n p (f l )£p(h)jp(h)tp(h)} of 
the nodes (n = 0, 1, . . .), orbital angular momentum, to- 
tal angular momentum, and isospin. The combined label 
{£j)a = ja + £ Q — 1/2 is used to label the £j combination 
uniquely; then £ = \{£j+l)/2] and j = [£j/2]+l/2 (where 
[x] is the integer part of x). The Greek letters a,(3, ... will 
be used to denote single particle states of either kind (p 
or h). 

The operator Q' x that creates an excited state |A) of 
energy E\ — Uu>\ with respect to the + ground state |0) 
and of angular momentum JM, 



\X;JM) = Ql JM \0), Qa,/a/|0) =0, 
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is written as 
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Henceforth the indices JM will be omitted, but implied 
throughout. 

The 2p2h state (and the corresponding creation oper- 
ator) is characterized, besides JM, by J p and Jh, the 
angular momenta to which the two particle states and 
the two hole states, respectively, are coupled. The same 
holds for the amplitudes X, y. Moreover, an ordering of 
the single-particle states is introduced and only the oper- 
ators with pi < P2, hi < h 2 are included in the expansion 
@ to avoid multiple counting of configurations. For ex- 
ample, the present convention is that a < (3 if t a < tp, or, 
for states of the same isospin, if (£j) a < {£j)p, or, when 
all other quantum numbers are the same, if n a <np. 

The SRPA ground state, which formally is the vacuum 
of the annihilation operators Q\, is approximated with 
the Hartree-Fock (HF) ground state. The latter is the 
Slater determinant that minimizes the expectation value 
of the given Hamiltonian, H = Hi+H 2 . The forward (X, 
X) and backward (Y, y) amplitudes are the solutions of 
the SRPA equations in ph © 2p2h space 
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The vanishing blocks are due to the choice of ground 
state. A and B are the usual N\ x N% RPA matrices 
(Ni the number of ph configurations), whose angular 
momentum-coupled forms are given by 



[^4-]p/l;p'fl' 



(e p — eh)S pp 'Shh' 

+ (ph- 1 ;J\H 2 \p'h'- 1 ;J) 

(e p — eh)S pp 'Shh' 

+ ^(-l)^+V--'i(2J 1 + l) 



Ji 



Plhip 2 h 2 J p Jh Plh1P2h2.Jp.Jh 
Pl<p 2 ,h 1 <h 2 ;J p ,J h 

\ I / j ■^Pih 1 p 2 h 2 Jp.J h Plhip 2 h 2 J p \}fJ 

Pl<p 2 ,h 1 <h 2 



l'l J ] \W-.J,\ii 1 \h l l-..h)m 



where O^ff^ creates a ph state and &p?hl P2 h 2 j j h creates 
a 2p2h state, coupled to the given quantum numbers. In 
particular, we have 

°pf = E i-^ jh - mh (jprn P j h -Tn h \JMj)al mp a hmh 

(3) 



[B] ph ., p , h , = ((ph-^J^p'ti K ,J)\H 2 \0) 

= E(" 1 ) J ' ,+v+J_ ' 7l ( 2Ji + 1 ) 

Ji 

x{f h % j}(l + V) 1/2 (l + ^') 1/2 

x(pp';J 1 \H 2 \hh';J 1 ). (7) 

The Ni x N 2 submatrix Ai 2 (N 2 the number of 2p2h 
configurations) describes the coupling between ph and 
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2p2h states, 

[•Al2]ph;p 1 p 2 h 1 h2.JpJh = 

(p/T 1 ; J\H 2 \( PlP2l J p )(h x h 2 ; A)- 1 ; J) 

= [i - {-iy h ^- Jh P{h u h 2 )]5 hl h 

x( _l)>+^ 2 +^+./ h (l + 5 hlh2 )- l ' 2 J p J h 
x{/ P 1 f)(PiP2;J p m P h 2 ;J p ) 

I Jli! Jh 2 Jp ) 

_ [i_ (-iy»i+^- J »p( PuP2 )}8 PlP 

x +^ +J+Jh {\ + Sp 1 p 2 )- 1 ' 2 JpJ h 

x{ i h ■ J A{hp2;J h \H 2 \h x h 2 ;J h ), (8) 

^ Jpi Jp2 Jh J 

while the N 2 x N 2 matrix A 22 contains the 2p2h states 
and their interactions, 

[■A 22 ]p 1 p 2 h 1 h 2 J P Jh;p'iP' 2 h 'i h '2- J 'p- J 'h = 
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xD-D^ + -(2i + i){J; #J 



propagation of a 2p2h state; the free propagation of the 
two particles while the two holes interact; the free prop- 
agation of the two holes while the two particles interact; 
and the free propagation of a ph pair while the other par- 
ticle and hole interact. In the following we will omit the 
J p , Jh indices for simplicity. 

If we neglect the coupling amongst the 2p2h states, A 22 
becomes diagonal and its elements are determined by the 
unperturbed 2p2h energies (diagonal approximation), 

[■A 22 ]p 1 h 1 p 2 h 2 ,p' 1 h' 1 p 2 h 2 — 

Sp 1 p' 1 S hlh ' i 5p 2 p' 2 S h2h ' 2 (e pl + e P2 - e hl - e }l2 ). (10) 

The validity of the diagonal approximation is examined 
in Sec. HVH 

As long as we are interested only in the single-particle 
response, determined by the ph amplitudes, X x and Y x , 
we may eliminate the X x , y x amplitudes from Eq. ([5]) 
and reduce the SRPA problem to an energy-dependent 
eigenvalue problem of the dimension of the RPA ma- 
trix 0], 



( A{E X ) B 
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In general, the expression for Aphp'h'(E) will involve the 
inverse of [(E + ir] 2 )X^ 2 — A 22 ], where 1n 2 the N 2 x N 2 
identity matrix. Within the diagonal approximation we 
have simply 



Aphp'h'{E) — Aphp'h' 

\A\ 2 \ph;p 1 p 2 hih2 \A\ 2 \p 1 p2hih2\p' 'h' 



E 
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x 



Jp Jp E 



Jh J'h E 



Jp[ Jpi Jp2 J l_ Jh[ Jhi Jh 2 

x^(-ir i+J » , r Ji (2j 1 + i) 

Jl 

x{j; j; ^(px^Jii^aiPiAuJi) 

In the above, e a are the HF single-particle energies and 
the operator P(a, (3) exchanges the indices a and j3. The 
two terms in A\ 2 describe the free propagation of the 
hole h while the particle p interacts with an intermediate 
ph state (terms with hi = h or h 2 = h, i.e., contain- 
ing 5h x 2 h), and the free propagation of a particle while 
the hole h interacts with an intermediate ph state (terms 
with S Pl 2P ). Thus self-energy corrections via bubble dia- 
grams are introduced either to a particle or a hole state, 
or a ph excitation is exchanged between a particle and 
a hole state. The four terms in A 22 describe: the free 



(9) 



, , E - (e pi + e P2 ~ e hl - e h2 ) + ir\ 2 

p 1 p 2 h 1 h 2 

A finite constant r/ 2 > is used in applications to 
smoothen the poles of the function and to introduce a 
width to the 2p2h states. 



A. Quantities of interest 

The quantities of interest are transition strength dis- 
tributions, or strength functions, Rp(E) of transition op- 
erators 

R F (E) = ]T|(A|Ft|0)| 2 5(£-£ A ) (13) 

A 

= ^5 F (£ A )<5(£-£ A ), (14) 



and their energy moments 



(15) 



determined, in general, by the amplitudes X , Y , X , and 
y through 



<A|Ft|0)=]T[/ p ,A p V-/p^ 



A * 
ph 



pit 
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Pipihih 2 
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where the coefficients /, / depend on the operator and 
ground state. Centroid energies can be defined as 



and 



E = mi/rriQ, 



(17) 



where the sums may be evaluated over the whole 
spectrum (unrestricted summations in Eq. (fT5|) ), or only 
in the energy region of a resonance. The width of a dis- 
tribution or a resonance can be expressed as 



A = 



m 2 
to 



E 2 . 



(18) 



Smoothened strength functions can be produced, for 
analysis purposes, by folding the discrete strength func- 
tions with a Lorentzian of width T, which yields 



Rf{E) 



r 



(£-£ A ) 2 + r 2 /4' 



(19) 



In the limit T — > 0, Eq. (|14p is recovered. 

We will consider IS and IV transitions of definite spin 
and parity J T , described by standard single-particle tran- 
sition operators Q. Then, for the HF ground state, we 
have 

Bf(E x ) = £[X p V + (-l) J Y p i*}{ P \\F\\h)\ 2 . 

ph 

(20) 

The two-body X, y amplitudes do not contribute. The 
energy moments mo and TOi will be the same in SRPA 
as in RPA ||. A (2 J + 1) multiplicity will be included 
in our final results. 

Two-body operators can also be considered. Then 
all amplitudes contribute to the transition matrix ele- 
ment. For example, let us consider the double dipole 
resonance, excited by the two-body operator Fddr.j = 
[F\yd ® -Fivd] j+ , where Fivd the usual (single-particle) 
isovector 1~ operator and J = 0, 2. We then have 



(A|4 DR;J |o) 



E[* P V + (-i) J ^UT R;J 

ph 



• / j [ PlP2hltl2 ~<~ 1 p 1 p 2 h 1 h2lJ piP2h! 

PiP2hih 2 



h 2 ' 



(21) 



where [To| 
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(23) 



Finally, we may define the strength distribution of a 
\pli~ 1 ) configuration, coupled to a given angular momen- 
tum state, via the quantity 



Sph(E) = J2(\ X ph\ 2 -\Yp\\ 2 )t(E-E x ) 
\>o 

= Y, s Ph(Ex)S(E~E x ), (24) 

A>0 

where the summation is over all eigenstates with E\ > 
0. (In the unperturbed case the centroid of S p h(E) is 
trivially identical to the HF ph energy e p — e/j and its 
width is zero.) Similarly, the strength distribution of a 
2p2h state can be defined using the X and y amplitudes. 
Energy moments and centroids, as well as smoothened 
distributions, can be defined as usual. 

We note that the total Too and toi (and centroid) of 
S p h will be the same in RPA and SRPA, since S p h is 
the sum of the strength functions related to the operator 

Oph 1 ^ an d its adjoint, too, in particular, should always 
amount to one. 



B. Related approximations and ground-state 
correlations 



(22) 



By setting the coupling matrices A12 and A22 and the 
2p2h amplitudes X, y equal to zero in Eq. ([5]), we re- 
cover the usual RPA problem. If, in addition, we ne- 
glect the ph residual interaction (i.e., B p h, P >h> = and 
A p h,p'h' = (e p - e h )S pp '5 h h'), we obtain a trivial, un- 
perturbed problem, where the eigenstates |A) are the ph 
configurations at the HF level and the Y amplitudes van- 
ish. 

By setting only B = in Eq. (J5j) , we obtain 
a second-order Tamm-Dancoff approximation (STDA), 
which amounts to solving the eigenvalue problem of the 
A block of the SRPA matrix. The backward amplitudes 
Y and y vanish and ground-state correlations implicitly 
taken care of by those are ignored. Setting also the cou- 
pling matrices A12 and A22 equal to zero, one gets the 
usual, first-order Tamm-Dancoff approximation (TDA). 

In a manner analogous to TDA, STDA is equivalent 
to a diagonalization of the Hamiltonian in the ph © 2p2h 
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space. It should be noted, though, that for J 71 " = + the 
HF ground state does not decouple from the STDA space 
(unlike TDA) . A diagonalization in the model space that 
includes in addition the HF state would produce a new 
ground state of lower energy. The + excitation spectrum 
would also be affected. 

In Sec. lIVEl we will investigate the possible influence of 
ignored ground-state correlations (GSC) on our results. 
A rigorous way to do that would be the use of an extended 
SRPA method built on a correlated ground state as self- 
consistently as possible (or other extended methods such 
as a self-consistent Green's function method [HI]). Since, 
however, that is a demanding project going beyond the 
scope of the present work, we resort instead to simpler 
approaches. One way to assess the role of GSC is to ig- 
nore them completely. This is accomplished within the 
(S)TDA. As a second approach, we have devised a sim- 
ple and rather heuristic renormalized version of SRPA 
(RSRPA), which takes into account to some extent the 
depletion of the Fermi sea. 

We start with a renormalized RPA approach (RRPA). 
Following the simplified RRPA method of Refs. [H Q]| 
- see also Ref. [14| — we assume partially occupied 
single-particle states and renormalize the residual cou- 
plings and the transition matrix elements accordingly. 
In particular, when calculating the matrix elements of 
A, the H2 terms in Eq. ^ are multiplied by a factor 



d [1] l =D 1/2 D 1/ , 2 , 

ph.p'h ph ph 



where 



D 



ph 



n h - n p 



(25) 



(26) 



and n a is the occupation probability of the orbital a. 
The matrix elements of B, Eq. J7J), are multiplied by the 
same factor. Finally, the single-particle transition matrix 
elements f p h, f p h, are renormalized by a factor D^ h 2 . 

In RSRPA we have to renormalize the matrix elements 
of A\2 and A22 as well. The former, Eq. ([8]), will be 
multiplied by the expression (f2"5)) . taking account of the 
occupation probabilities of the p, h states, as well as by 
a factor 

[2] 1 

to take into account the occupation probabilities of 
the pi,2, ^1,2 states. Similarly, the H2 terms in the 
A22 matrix elements, Eq. will be multiplied by 



,"V.' h' h' ■ 



Pih 1 p 2 h 2 p^h^p'^h'^ 

In practice, we will not solve these equations iter- 
atively. The single-particle energies and eigenstates 
needed to evaluate the above matrix elements will be the 
HF ones, whereas the occupation probabilities are calcu- 
lated using the shell model. Obviously, this method is 
neither rigorous nor consistent, but it should help us get 
an idea regarding the sensitivity of our results to GSC. 



C. Consistency and stability 



It is well known, that in self-consistent RPA (meaning 
that the same Hamiltonian is used to calculate the HF 
ground state and the residual interaction) the spurious 

state related to the CM momentum operator P (single 
spurious state) will appear at zero energy and be exactly 
separated from the physical spectrum, provided that all 
ph and hp configurations available in the (sufficiently 
large) single-particle space are taken into account. The 
same does not hold, however, when extensions of RPA 
are considered, whether that means considering a corre- 
lated ground state or higher-order configurations. As was 
formally shown in Ref. [l5(], once 2p2h configurations are 
included in the model space, a necessary condition for 
the single spurious state to appear at zero energy is that 
all single-particle amplitudes be taken into account: This 
means not only ph and hp amplitudes, but also pp and 
hh. It was in fact shown that pp and hh amplitudes affect 
the spurious state even with a HF ground state, because 
at energies equal to e p — e p > or — e/j/ they do not van- 
ish. Obviously our approach does not include them. This 
problem will be examined quantitatively in Sec. IIVEI 



A more severe problem in HF-based SRPA is the onset 
of instabilities. Contrary to RPA, the "self-consistent" 
use of the HF ground state does not guarantee that the 
SRPA stability matrix will be positive-definite: Thou- 
less's theorem was proven specifically for ph excita- 
tions [16], [l7| . Thus RPA is related to the stability condi- 
tions for the HF solution, but SRPA seems related to an 
extended variational problem @, EH. In practice, we will 
find that low-lying states appear at imaginary or nega- 
tive energies. 1 It is not difficult to demostrate how neg- 
ative eigenvalues can occur, already at the STDA level 
(B = 0), if we consider Eq. (fT2]) : For energies E below 
typical 2p2h energies, strong (or many) A12 elements, 
regardless of sign, can cause even the diagonal elements 
of A(E) to become negative and destroy the positive- 
defmitness of the matrix. Additional ground-state cor- 
relations could cure this problem, through a renormal- 
ization of the A12 couplings, as illustrated above, or by 
filling the vanishing B— submatrices of Eq. |5]| with finite 
elements (see also Ref. [ll[ for a related discussion within 
an extended dressed RPA). 



The question whether SRPA gives meaningful results 
for giant resonances despite the problematic solutions at 
low energies is also tackled in Sec. IIV El 



1 Here we speak of a "negative-energy" excitation when the eigen- 
state with E\ < is normalized to +1, i.e., the norm of the 
positive-energy counterpart is negative. 
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III. SOLVING THE SRPA PROBLEM 

In practice we proceed as follows. We first choose 
a single-particle space, consisting of harmonic-oscillator 
eigenstates. The larger the space, the better the conver- 
gence of our results to their final values. All angular- 
momentum coupled two-body matrix elements of the 
given interaction within the harmonic-oscillator basis are 
calculated in advance as described in Ref. Q and stored. 
We solve the HF equations within the given space and 
then we use all ph and 2p2h configurations, which are 
available within the space and can couple to the desired 
angular momentum and parity, to construct the SRPA 
matrix. The single-particle space is characterized by ei- 
ther 

• £max = (2n + f) m „, the number of energy quanta 
in the highest oscillator shell considered. This 
means that all single-particle states within the low- 
est e max + 1 shells are used. An additional cutoff 
f max in £ may also be imposed. 

or 

• rimax and £ max . A cutoff £ max in £ is considered 
and each possible £j state is expanded in the lowest 
"max + 1 harmonic-oscillator £ states. 

Since angular momentum is preserved, for an RPA cal- 
culation of the J 71 response, one has to include j's up to 
the maximal occupied one plus J, so that all j p 's with 
\jh — jp\ < J < jp + jh are included in the space, and 
may omit higher j p 's. In SRPA there is no j-cutoff pro- 
vided by angular momentum conservation. In principle, 
particle states with infinitely large j's can contribute to 
allowed 2p2h configurations and such configurations can 
couple with ph ones through finite matrix elements of 
Ai2- In particular, infinitely large j Pl and j P2 can cou- 
ple to any given J p , thus resulting in finite values of 
(P1P2; J P \H 2 \ph 2 ; J p ) — see r.h.s. of Eq. (JH). For re- 
stricted n max , however, upper bounds are provided by 
the transformation (Moshinsky) brackets needed to eval- 
uate the above matrix element in the HO basis [19]. Nat- 
ural cutoffs on n max or e max should be provided by the 
properties of the interaction. Perhaps such considera- 
tions can be used, along with the convergence behavior 
of the SRPA solutions with respect to the single-particle 
basis, for the optimization of the latter, but so far this 
has not been done. 

The next task is to solve the SRPA eigenvalue prob- 
lem. First we note that all the submatrices comprising 
the SRPA matrix, Eq. (O, are real in the cases studied 
here. The submatrices A, B, and A22 are also symmetric. 
The same holds for the N x N blocks indicated in Eq. ([5]) 
(separated by lines). The dimension 2N of the SRPA ma- 
trix is given by twice the number of ph and of 2p2h con- 
figurations available in the model space (JV = N\ + N 2 ). 
The latter, N 2 , can be rather large: For the purposes of 
the present work we encountered problems with N up to 
10 6 . In the following we present technical information 



on how we deal with the SRPA eigenvalue problem in 
practice. We shall distinguish three cases: 

• All eigenvalues and eigenvectors are calculated. 
This is feasible for relatively small problems, e.g., 
N ~ I0 4 . 

• For large problems the complete solution becomes 
impossible. Then only a small portion of consecu- 
tive eigenvalues and the respective eigenvectors are 
calculated. In most cases, all excitations of interest 
lie at the lower end of the spectrum, so it suffices to 
evaluate only the lowest (a few tens or hundreds) 
positive eigenvalues. 

• No eigenvalue problem is solved. The response 
function is evaluated directly. The strategy is well 
known and very practical when one is interested in 
the single-particle response, and only in the final 
result, namely a smoothened strength function. 

1. Small problems: all eigenvalues 

Standard library routines can be used to solve small 
enough eigenvalue problems. The solution can be sped 
up considerably if the SRPA problem is reduced to half 
the dimension, N x N. A method relying on a Cholesky 
decomposition of the matrix A + B or A — B is the most 
common way to do it, with the additional bonus that it 
produces a symmetric N x N eigenvalue problem [20| . 
However, the method will not work if neither A + B nor 
A — B is positive-definite, a problem that can occur when 
there are instabilities, or just spurious states at imaginary 
energies. We employ a modified method instead, relying 
on a generalized Cholesky decomposition [2l|, which is 
equivalent to the original one when the matrix under de- 
composition is positive definite (thus still producing a 
symmetric eigenvalue problem) and still works when it is 
not (producing a non-symmetric eigenvalue problem of 
half the dimension), with minimal additional computa- 
tional effort. For details, see Ref. [2l| . 

2. Large problems 

When the dimension of the matrix is very large, the 
first problem that occurs is the storage of the matrix el- 
ements. Fortunately, most of the SRPA matrix matrix 
elements vanish, especially when the diagonal approxi- 
mation is used, Eq. (fTU|) , so it becomes possible to store 
all the finite ones in memory. Obviously, one needs to 
store only one N x N matrix (e.g., A or A ± B) and one 
small, Ni x Ni matrix (B). Moreover, only the upper 
(or lower) triangles need be stored, since the matrices 
are symmetric. 

The dimension and degree of sparseness of the SRPA 
matrix depends on the nucleus, type of response, and 
model space. As an example, let us mention that in the 
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case of the 2 + response of 40 Ca with the Vjjcom, and for 
spaces large enough for reasonable convergence, there are 
about a million configurations. 0.5 — 1GB of storage are 
needed in the diagonal approximation, but 100 GB could 
be required for the full problem (10 10 elements in double 
precision) . 

Again the dimension of the SRPA problem can be re- 
duced by half, to speed up the numerical solution. In 
order to save matrix operations we chose not to perform 
a decomposition, but resort to a more straightforward re- 
duction method [2l| : We solve the non-symmetric N x N 
problem 

(A -B)(A + B)R X = E 2 X R X , (28) 

where the eigenvectors R x = E x 1 ^ 2 (X X + Y x ) obey the 
normalization condition (R X ) T (A + B)R^ = ±<5a m (real 
and positive E\). The properly normalized X and Y 



arrays (|A A | 2 — 



\Y X 
1 



= ±1) are then given by 
1 



Y A = - 




(29) 
(30) 



where In is the N x N identity matrix. 

Finally, an Arnoldi iteration procedure from the 
ARPACK package [22j is employed to solve the prob- 
lem (|28p only for the k lowest positive eigenvalues E x , 
where k << N. In principle, it is possible to solve for 
the first k eigenvalues lying above a given energy value 
E 2 Ssct , not necessarily equal to zero. In practice, how- 
ever, such a strategy can prove problematic: Already at 
moderate energies the density of eigenstates can be so 
large (see Sec. IIV[) . that the Arnoldi iteration will likely 
fail to converge. 



3. Response function formalism 

As long as one-body transition operators are consid- 
ered, we may just solve the reduced SRPA problem, 
Eq. (JTTJ) . In general, the reduction procedure involves 
the inversion of a large matrix (of the dimension of the 
2p2h space), but when A22 is diagonal, that is reduced 
to a trivial number inversion, see Eq. (|12p . There are 
ways to solve such an energy-dependent eigenvalue prob- 
lem [23|, [24l ]. An efficient alternative is to employ the 
response-function formalism. Then, instead of explicitly 
solving the eigenvalue problem, one can obtain directly 
the strength function of interest 0, HH]. First, the ph 
Green's function, a iVi x N\ matrix, is evaluated, 



G(E) 



where the elements of the Aq — dimensional (assumed 
real) array T (JF), in ph space, are the matrix elements 
of the transition operator, f p h (—fhp) — see Eq. (fTBT) . 
Finally, the strength function of interest is given by 



A{E) - (E + irn)l Nl B 

B* A*(-E) + {E-i m )I Nl 

(31) 

(Vi ~~ * + ) by numerical matrix inversion. The response 
function for a given single-particle field is given by 

U(E) = (T T , T T )G(E) ( J V (32) 



R(E) 



: 11{E) 



(33) 



In practice, the Green's function, response function, 
and finally the strength function are evaluated over the 
energy range of interest, which is represented by mesh 
points Ei. The constants 771, which provide the ph and 
2p2h states, respectively, with a finite width should, 
in principle, be small enough for the structure of the 
strength function to be resolved as desired. The mesh 
size SE = Ei+i — Ei should be smaller than both of 
them (a factor 3-4 suffices) . The choice 771 = 772 produces 
a smoothened strength function that is practically the 
same as the discretized strength function (obtained by 
explicit diagonalization of the SRPA matrix) folded with 
a Lorentzian of width T — 2rji 2, cf. Eq. (flTJl) . 



IV. RESULTS 

In the following, we will discuss the features of the 
SRPA solutions with the help of illustrative examples. 
We note that we will use the acronym SRPAO when refer- 
ring specifically to the diagonal approximation, Eq. (|10[) . 
and SRPA(O) when referring explicitly to both solutions, 
with and without the diagonal approximation. The in- 
trinsic nuclear Hamiltonian is employed. In all cases we 
set the oscillator length parameter of the single-particle 
basis equal to b = 1.7 fm (for Vucom) or b = 1.8 fm (for 
V BB ). 

It will prove instructive to consider, among others, a 
relatively small SRPA problem, for which we can cal- 
culate all eigenstates — with and without the diagonal 
approximation. The purpose of such an exercise is not 
to perform a realistic calculation, but to illustrate and 
discuss selected features of the method: in particular, 
how the large amount of 2p2h configurations influences 
the response function and the distribution of the SRPA 
eigenstates. 

As such a "toy model" we choose the monopole (0 + ) 
response of 16 in a rather small single-particle space, 
consisting of seven oscillator shells, and using the Vucom 
interaction. In total, there are 4148 positive-energy 
Eigenstates to evaluate (Aq = 14, ^2 = 4134). The 
single-particle, isoscalar monopole strength distribution 
is shown in Fig. [1] in linear (Fig. QJi) and logarithmic 
(Fig. [Tb) scale. 

Further examples will be introduced in the following 
as needed. 
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FIG. 1: Isoscalar monopole response of le O, calculated within a single-particle basis of 7 oscillator shells: Strength function 
calculated within SRPA and SRPAO (diagonal approximation, Eq. (HSJ), a s well as RPA, in (a) linear and (b) logarithmic scale. 



A. Diagonal approximation 

When looking at Fig. [TJi, we observe that the results 
of the diagonal approximation, SRPAO, are very close to 
the exact SRPA results. It looks as though the effect 
of the 2p2h space on the single particle response is ap- 
proximately the same, whether or not the 2p2h states 
are considered unperturbed. We have verified that the 
approximation is equally good in larger spaces and for 
different types of response — always of single-particle 
operators. See, e.g., the IVD response of 16 calculated 
within a space of 13 oscillator shells, as shown in Fig. 1 
of Ref. The approximation remains quite good for 
heavier nuclei, for example 48 Ca, whose ISQ response is 
shown in Fig. [2^, calculated within a space of 9 shells. 

We expect the diagonal approximation to be reliable 
for soft, perturbative interactions in general. The addi- 
tional couplings within the 2p2h space, ignored in the 
diagonal approximation, constitute higher-order correc- 
tions to the excitation propagator with respect to the in- 
teraction, as can be demonstrated diagrammatically [2|| . 
Results with the Brink-Boeker potential, Vbb, corrobo- 
rate this speculation. Vbb is even softer than the Vucom- 
For example, it produces much smaller second-order cor- 
rections to the nuclear binding energies in the perturba- 
tion expansion beyond Hartree-Fock [26] . In Fig. ^> we 
show the IVD response of 16 O calculated using Vbb- The 
diagonal approximation is very good in this case. 

To be on the safe side, it is always advisable to verify 
the quality of the approximation within some tractable 
space. In Ref. [27], for example, it was found that it is 
quite bad in the case of metallic clusters when the bare 
Coulomb interaction is used. 

From Fig. QJi we deduced that the SRPA and SRPAO 
solutions yield almost the same results for the giant res- 
onance. However, when we look at the strength distribu- 
tion on a logarithmic scale, Fig. [TJa, it becomes obvious 
that the two solutions give different results for the ma- 
jority of eigenstates, which are mostly of 2p2h character. 
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FIG. 2: Quality of the diagonal approximation, (a): IS 
quadrupole response of 48 Ca in a single-particle basis with 
tmax = 8 and using Vucom- (b): IV dipole response of 16 O 
in a single-particle basis with e max = 12 and Imax = 8 and 
using Vbb- Results are shown obtained with RPA, solving 
the full SRPA problem and using the diagonal approximation 
(SRPAO). (In all cases, P = 0.5 MeV.) 



One implication is that the diagonal approximation can- 
not be relied upon when examining, e.g., double giant 
resonances and the response function of two-particle op- 
erators in general. As an example, in Fig. [3] we plot 
the strength function of the + component of the double 
dipole resonance. We notice that the SRPAO strength 
function is very close to the unperturbed one (HF) and 
differs significantly from the SRPA result. 

In Fig. 5] we show the density of eigenstates (number 
of states per 5 MeV excitation energy) when solving the 
SRPA problem with and without the diagonal approx- 
imation, as well as the density of unperturbed ph and 
2p2h states (HF). One immediately notices that for the 
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FIG. 3: + component of the double dipole resonance of le O, 
calculated within a single-particle basis of 7 oscillator shells: 
Strength function (r = 0.5 MeV) calculated within SRPA and 
SRPA0 (diagonal approximation, Eq. (|10|) '). as well as using 
the unperturbed (HF) ph and 2p2h states. 
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FIG. 4: For the + response of ls O, calculated within a single- 
particle basis of 7 oscillator shells: Number of SRPA and 
of SRPA0 (diagonal approximation, Eq. (|10|) eigenstates per 
5 MeV excitation energy — in logarithmic scale. The corre- 
sponding density of unperturbed (HF) ph and 2p2h states is 
also shown. 



most part the SRPA0 density of states is very similar to 
the unperturbed one. 

Finally, below about 40 MeV, where the giant reso- 
nance lies, the SRPA and SRPA0 densities practically 
coincide. This is not surprising, given the very small 
number of eigenstates present there, but it is not always 
the case. We note, for example, that in the case of the 
2 + response of 48 Ca, where there are several unperturbed 
2p2h states at low energies, the SRPA0 density of states 
follows mostly the HF one in that energy region. The 
single-particle strength functions, however, are almost 
identical (Fig. [2^,) and the diagonal approximation re- 
mains well justified. 



B. Downward shift of resonances 

In Figs. [T] and [5] we notice that the GR lies lower in 
SRPA(0) than in RPA. This is a rather general result — 



see also Ref. Q . It is a general feature of the coupling 
to 2p2h excitations [25|, HH (and not particular to nu- 
clei [27|), and shows why traditional effective interactions 
cannot be used in SRPA. The self-energy corrections re- 
sponsible for the modification of the ph energies and the 
lowering of the resonance energies (see also Sec. lIV Cj) are 
already parameterized in the interaction, so that realistic 
results are obtained already at the RPA level. Employing 
such an interaction in SRPA will result in double count- 
ing of those effects. 

Ths problem is circumvented in practical applications 
in the literature through the use of realistic single-particle 
energies and subtracting procedures, which remove the 
real part of the self energy, responsible for the shift. 



C. Fragmentation in SPRA 

As a microscopic theory of collisional damping, SRPA 
has been used extensively to describe the spreading width 
and strength fragmentation and quenching of collective 
excitations. Let us examine our results in this context. 

No spreading width is observed in Fig. as there 
are no configurations available in the vicinity of the res- 
onance. But when we look at Fig. [TH, we realize that 
there is practically a continuum of weak SRPA or SRPA0 
eigenstates (N = 4148) extending to high energies. Most 
are predominantly of 2p2h nature and their contribution 
to the single-particle strength is rather small. Neverthe- 
less, they carry a non-negligible percentage of the total 
strength and provide a mechanism of strength quench- 
ing for the giant resonance (by contrast, there are only 
N% = 14 RPA eigenstates). We note, in particular, that: 
(i) As expected, the total energy-weighted strength mi is 
found to be practically the same in all three cases (RPA, 
SRPA0, SRPA), and the total strength mo is found just 
about 3% larger in SRPA or SRPA0 than in RPA, but (ii) 
the strongest ISM peak appears at 21.36 MeV in the case 
of RPA, 16.50 MeV in SRPA0, and 16.25 MeV in SRPA 
and (iii) the six lowest eigenstates carry a total strength 
of 313.4 fm 4 in RPA, which is almost all the RPA strength 
(they lie at 21 to 47 MeV excitation energy), 304.5 fm 4 
in SRPA0 (16 to 34 MeV) and 304.0 fm 4 in SRPA (16 to 
34 MeV), i.e., about 3% less in SRPA(0) than in RPA; 
finally, (iv) the width of the ISM strength functions be- 
low 40 MeV is 3.86 MeV in RPA and 3.73 MeV in SRPA, 
while the widths over the whole spectrum are 4.65 and 
18.94 MeV respectively. 

The coupling with 2p2h configurations affects not only 
the collective excitations, but at the same time the single- 
particle states. Within the present formalism one can 
demonstrate this through the strength of the ph configu- 
rations. In Fig. O for example, we plot the HF, RPA 
and SRPA strength distribution S ph (E), Eq. ([23]). of 
the ph configurations \{yPz/2){ v ^P3,/2) 1 ] + )> contribut- 
ing to the + strength of Fig. [TJ In HF, these are 
well-defined transitions at energies equal to e p — eh with 
strength 1. In RPA the ph transitions appear fragmented. 
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FIG. 5: Fragmentation and shift of ph states — + response 
of 16 in a space of 7 shells. Thin dark bars show how the 
spectroscopic strength S p h{E), Eq. (|24|l . of the ph configura- 
tions 1(^3/2X^(^3/2) + ) (contributing to the monopoole 
strength) is distributed within (a) SRPA, (b) RPA and (c) 
HF. Thicker, pale (cyan) bars denote the distribution of 
|(^lp3/ 2 )(^0p3/2) _1 ; + ) (one particle shell only). 



There is a cluster of configurations at energies around 
30 MeV, which can be attributed to the 2hu> transition 
|(^lp3/2)(^0p 3 / 2 ) _1 ; + ). A second cluster of states is 
visible around 60 MeV roughly corresponding to Ahu) con- 
figurations, and so on. In SRPA the strength distribu- 
tion appears even more fragmented, as well as shifted to 
lower energies. The shift reflects an effective compression 
of the single-particle spectrum (a modification of the nu- 
cleon effective mass), with respect to HF, and leads to the 
downward shift of the collective states, discussed above. 

Although we can visually identify the different shells 
contributing to the strength distribution, and observe an 
energetic shift, the total mo and m\ — and thus the cen- 
troid — of any given ph configuration is the same in RPA 
and SRPA. Nevertheless, in RPA we find only 1% of the 
|(^lp3/2)(^0p 3 / 2 ) _1 ; 0+) m strength above 40 MeV. In 
SRPA the strength in the same region is 14%. When 
we look at the distributions in logarithmic scale, Fig. [51 
we realize that the ph strength distributions span the 
whole space available, reflecting the extended spectro- 
scopic functions of the particle and hole states. 



D. Truncation procedures 

As is obvious from Figs. Q] and [4] the majority of basis 
states (unperturbed ph and 2p2h states) lies at high ener- 
gies relative to the giant resonance. Figure[7J(see caption 
for details) illustrates what would happen if we truncated 
the model space by setting an upper cutoff I?2p2/i,max to 
the 2p2h energies taken into account. Had we set, e.g., 
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FIG. 7: For the + response of ls O, calculated within a single- 
particle basis of 7 oscillator shells, but imposing an energy 
cutoff on the 2p2h configurations taken into account, Sph.maso 
Energies of the four lowest eigenstates (axis on the left) and 
percentage of basis states taken into account (axis on the 
right), vs -Eph.max- The arrow indicates the point at which 
all available states are used (maximal 2p2h energy). 



-E2 P 2/i,max = 100 MeV, we would have excluded 75% of 
the configurations available in this particular space and 
the eigenenergies would have been at least 2 MeV higher 
than their converged (with respect to -E^h.max) values. 
Inclusion of the energetically lowest 50% of configurations 
would still not warrant good convergence. 

We found that the convergence of the eigenenergies and 
the saturation of the model space follow the same pattern 
as demonstrated in Fig. [7J regardless of model-space size 
and type of response. 

Finally, we should note that most finite A12 elements 
are small. The distribution of their values is always cusp 
shaped around zero and a vast majority of them have am- 
plitutudes no larger than 5% the value of the strongest 
element. It would be computationally economical to ne- 
glect those elements. We found, however, that such a 
procedure influences the position of the resonances as 
well as the shape of the strength distribution and cannot 
be blindly trusted. 
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TABLE I: Isoscalar 2 + spectrum of 48 Ca calculated in a basis 
of 9 oscillator shells: Energy E (in MeV) and strength B (in 
fm 4 ) of the two main low- lying states and of the giant reso- 
nance peak (GQR) within RPA and SRPA, "renormalized" 
RPA and SRPA (RRPA, RSRPA — see text), and TDA and 
STDA (SRPA with B = 0). (Whenever the GQR is split into 
two or three major peaks, the centroid and total strength is 
given.) 



E. Consistency and stability issues 

The self-consistent RPA that we have used produces a 
spurious state at practically zero energy, and leaves the 
rest of the spectrum uncontaminated [8j], as expected. 
The SRPA dipole spectrum, however, may contain spu- 
rious admixture, as discussed in Sec. Ill CI In order to 
quantify this problem, we have examined the IS dipole 
response. We found that relatively strong spurious states 
appear mostly at about 5 to 8 MeV. Using a transition 
operator of the usual radial form (oc r 3 — §(r 2 )r) and its 
uncorrected form (oc r 3 ), we found that only the lowest 
part of the dipole spectrum is strongly affected by the 
choice of operator, while there are no significant contam- 
inations in the spectrum around and beyond the IV GDR 
peak. An example is shown in Fig. [5] 

Let us note, that including a small percentage of 2p2h 
states would not shift the spurious state far away from 
zero, but for large spaces such as those used here the ef- 
fect is noticeable. As we couple the ph states to more and 
more 2p2h configurations, the spurious RPA state moves 
away from zero and may occur at imaginary, or even 
negative energies (or, equivalently, positive energies, but 
assuming negative norm — this is the case in Fig. [5]). It 
may also appear fragmented, as spurious admixtures en- 
ter the rest of the spectrum, or as additional eigenstates 
occur in its proximity. 

As anticipated in Sec. Ill CI another problem with the 
HF-bascd SRPA is the onset of instabilities. We find, in 
particular, that low-lying states (OTilu 2 + states and the 
collective 3~ excitation) appear at imaginary or negative 
energies. This is the case, for example, in Fig. (5^: the 
norm of the collective quadrupole state at 4.4 MeV is, in 
fact, negative; the positive-norm counterpart appears at 
—4.4 MeV. As discussed previously, the problem seems to 
be the inadequate treatment of ground-state correlations. 

Obviously, SRPA is not appropriate for describing low- 
lying states. The question is then whether it still gives 
meaningful results for giant resonances. We have tested 



the sensitivity of our results to CSC with the help of the 
renormalized SRPA devised in Sec. Ill Bl and of STDA 
(SRPA setting B = 0) and, in first order, the correspond- 
ing RRPA and TDA. Examples are shown in Fig.[9]for the 
octupole response of 16 O and Table U for the quadrupole 
response of 48 Ca. In all cases the diagonal approximation 
has been employed. 

The general trend can be described as follows: We be- 
gin with minimal CSC, in STDA. Inclusion of the B ma- 
trix and the backward amplitudes (SRPA) pushes the 
solutions to somewhat lower energies and in general af- 
fects their strength. The downward shift, which can lead 
to imaginary solutions, was interpreted recently in the 
framework of random-matrix theory [29j : The matrix 
B couples the positive-and negative-energy branches of 
the (S)RPA solutions and causes an attraction between 
them. Strong enough coupling leads to a merging of the 
two branches at zero energy and eventually to imaginary 
or complex solutions. Finally, we find that subsequent 
renormalization of the matrix elements (RSRPA) shifts 
the solutions back to higher energies, i.e., slightly closer 
to the RPA solutions, plausibly because of the weaker 
A12 couplings. 

The effect of renormalization appears rather small in 
the case of giant resonances and higher-lying solutions 
in general, but this is not the case for the low- lying 
quadrupole and octupole states: their energy shifts in 
RSRPA, relative to RPA, are noticably more moderate 
than in SRPA. The STDA results confirm in general that 
the lower-lying states are more sensitive to CSC than 
the higher ones. We note that, since STDA is a Her- 
mitean problem, it will eliminate the imaginary solutions 
of SRPA by construction — see, e.g., Fig. [5k. As we ob- 
serve in Table HI however, the same does not hold for the 
negative solutions. It is thus confirmed that the strong 
resonance shifts (with respect to RPA) in this case are 
induced by the A12 couplings, as discussed in Sec. Ill CI 

In most cases the RSRPA solutions lie higher than 
the STDA, i.e., renormalization has a stronger effect on 
the energies, with the notable exception of the giant 
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FIG. 9: Exploring the influence of ground state correlations: Isoscalar 3 response of 16 O (e max = 14, ^ max = 10) within RPA 
or SRPA and "renormalized" RPA or SRPA, (RRPA, RSRPA — see text), as well as TDA and STDA (SRPA with B = 0). 
(a): Evolution of the three lowest (closest to zero) eigenvalues as we increase the number of 2p2h configurations considered, 
starting from zero and up to all 2p2h states available in the single-particle space, N2. The numbers assigned to the various 
points indicate the corresponding 2p2h— energy cutoff imposed, i?2 P 2h,max, in MeV. The energy of the lowest eigenstate becomes 
imaginary for larger N2 in (R)SRPA; then its amplitude is indicated. Note that the STDA and SRPA results for the second 
state almost coincide, (b): Strength function for £ 2p 2h,max = 150 MeV in SRPA, RSRPA, and STDA (r = 2 MeV). (c): 
Strength function in RPA, RRPA, and TDA. 



monopole resonance (not shown), where the backward 
amplitudes seem more relevant. 

Finally, regarding our first-order results, renormaliza- 
tion affects the RPA results very weakly. The same con- 
clusion was reached in Ref . [3] > where a more consistent 
RRPA was applied. The backward amplitudes, though, 
missing in TDA, are found important for the description 
of low-lying states. 

The present results suggest that giant resonances (at 
least for J > 0) in SRPA are only moderately sensi- 
tive to the treatment of GSC. Of course, a comprehen- 
sive inclusion of GSC in the SRPA formalism would in- 
volve filling up the B— sections of the SRPA matrix with 
many more finite elements. It cannot be predicted how 
strongly those could affect the results. In the applications 
shown in Ref. (l8j |. the so-called extended SRPA includ- 
ing correlations did not produce strong corrections to the 
SRPA strength functions, especially after renormalizing 
the ground state [13, HH . 



V. CONCLUSIONS 

The present work constitutes a feasibility and justifi- 
cation study of large-scale SRPA calculations. The mo- 
tivation was the prospect of studying nuclear collective 
excitations using SRPA and unitarily transformed nu- 
clear Hamiltonians, which have not been fitted to first- 
order RPA results. The discussion, however, has been 
kept general. We showed how the large model spaces of 
SRPA can be treated and discussed the salient features of 
the solutions, including the energetic shift relative to the 
RPA solutions and the fragmentation of strength. The 



diagonal approximation was found reliable for soft inter- 
actions and for single-particle strength distributions. 

We found that low-lying states become unstable in 
SRPA, due to an inadequate treatment of ground-state 
correlations. Nevertheless, giant resonances and higher- 
lying solutions in general do not appear sensitive to 
ground-state correlations. We have thus concluded that 
SRPA can be applied in the giant-resonance region with 
reasonable confidence. 

SRPA is primarily the theory of collisional damping. 
A more comprehensive method to study nuclear collec- 
tive states should consider also coupling to collective low- 
lying phonons [12, [HI • It is not easy to tell at this point 
how strongly those could enhance or (partly) cancel the 
effect of coupling to the large amount of 2p2h states con- 
sidered in SRPA. In Ref. [28| it was shown that the res- 
onance shifts due to the 2p2h states are stronger than 
those due to collective phonons, although both mecha- 
nisms are expected to contribute to the damping width. 

Notwithstanding some particular shortcomings of the 
present SRPA formalism, this work paves the way for sys- 
tematic studies of giant resonances using finite-range in- 
teractions and extended RPA theories without adjustable 
parameters and arbitrary truncations of the model space. 
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